Thermal fluctuations and phase diagrams of the phase field crystal model with pinning 
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We study the influence of thermal fluctuations in the phase diagram of a recently introduced 
two-dimensional phase field crystal model with an external pinning potential. The model provides 
a continuum description of pinned lattice systems allowing for both elastic deformations and topo- 
logical defects. We introduce a non-conserved version of the model and determine the ground-state 
phase diagram as a function of lattice mismatch and strength of the pinning potential. Monte Carlo 
simulations are used to determine the phase diagram as a function of temperature near commen- 
surate phases. The results show a rich phase diagram with commensurate, incommensurate and 
liquid-like phases with a topology strongly dependent on the type of ordered structure. A finite-size 
scaling analysis of the melting transition for the c(2 x 2) commensurate phase shows that the thermal 
correlation length exponent v and specific heat behavior are consistent with the fsing universality 
class as expected from analytical arguments. 
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I. INTRODUCTION 



dard molecular dynamics^ 



Many two-dimensional (2D) lattice systems on periodic 
potentials form different ordered structures at low tem- 
peratures which may be unstable against thermal fluc- 
tuations and the mismatch between the competing peri- 
odicities. Adsorbed layers on crystal surfacesii 2 ^, vortex 
lattices in 2D superconductors^ and colloidal crystals^ 
in a periodic substrate are important examples of current 
interest. In numerical simulations, pure elastic models or 
particle models are often employed. However, there are 
important limitations in these approaches. In the for- 
mer, plastic deformations such as dislocations are not 
taken into account and in the latter, the accessible time 
and spatial scales may be very limited. These limita- 
tions are particularly important in atomistic simulations 
of crystalline materials described by microscopic models 
with complicated interactions, where the time scale is set 
by the vibrational frequency and consequently relatively 
small system sizes are often used in numerical simula- 
tions. 

Recently, a phase field crystal (PFC) model for pinned 
lattice systems was introduced^ that allows for both elas- 
tic and plastic deformations of the lattice. The model 
describes the lattice system as a continuous density field 
in an external pinning potential. In this formulation a 
free energy functional is introduced which depends on 
the phase field vb(r) corresponding to the particle den- 
sity. In absence of a pinning potential, the free energy is 
minimized when vb is spatially periodic? forming a hexag- 
onal solid phase in 2D. By incorporating phenomena on 
short length scales the model naturally includes elastic 
and plastic deformations and it should be computation- 
ally more efficient than particle simulations such as stan- 



In a previous worki, the phase diagram of the PFC 
model in a periodic pinning potential with square symme- 
try was determined as a function of the pinning strength 
and lattice mismatch, in the absence of thermal fluctu- 
ations. A numerical minimization procedure was used 
to determine the different ordered phases and the phase 
transitions, based on the corresponding dynamical equa- 
tion of motion for a conserved field. In the present work, 
we consider the influence of thermal fluctuations and lat- 
tice mismatch in the PFC model with pinning. We first 
introduce a non-conserved version of the model which 
is more efficient for equilibrium simulations, such that 
a constant chemical potential is included in the free en- 
ergy to fix the average density. With this new version, 
we determine the ground-state phase diagram as a func- 
tion of lattice mismatch and strength of pinning potential 
using the corresponding non-conserved dynamical equa- 
tions of motion. Monte Carlo (MC) simulations are then 
used to determine the phase diagram as a function of 
temperature and mismatch near commensurate phases. 
The results reveal that the model has a rich phase di- 
agram, which includes commensurate, incommensurate 
and liquid-like phases with a topology dependent on the 
type of ordered structure. In particular, a finite-size scal- 
ing analysis of the melting transition for the c(2 x 2) 
commensurate phase shows that the thermal correlation 
length exponent v and specific heat behavior are consis- 
tent with the Ising universality class, which is expected 
from analytical arguments based on symmetry consider- 
ations and Landau free-energy expansion. 
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II. MODEL 

The mean field free-energy functional of the PFC 
model with an external pinning potential^ can be written 
in dimensionless units as 

F = F J dx{^[r + (1 + V 2 ) 2 ]V + ^ + Vtf)}, (1) 

where ip(x) is a continuous field representing the local 
number density of the particles and V(x) represents the 
external pinning potential. The overall constant Fq sets 
the energy scale and r is a parameter. In Eq. ([1]), ifj{x) 
is a conserved field and the minima of the free energy 
depends on the average value of ip(x). To remove this 
conservation constraint on ip(x), we introduce here an 
additional linear term to the free energy containing the 
chemical potential //, which controls the average value 
of ip{x). This leads to a model mean field free-energy 
functional 

F = F J dx{^j[r+(l + V 2 ) 2 ]^ + ^ + V^-^}- (2) 

This modification allows the use of a MC algorithm with 
non-conserved dynamics as described in Sec. 3, which is 
more efficient for equilibrium simulations. 

In the absence of the pinning potential, V(x) = 0, the 
free energy of Eq. ([2]) can be minimized by a configu- 
ration of the field ip(x) forming a hexagonal pattern of 
peaks with wavector of magnitude fco — 1 , when the val- 
ues of the parameters r and \x are chosen appropriately. 
This structure of peaks can be regarded as a crystalline 
system and the free energy functional of Eq. (J2J) can then 
be used to describe both elastic and plastic properties of 
such a lattice system within a mean field description^^ 
for temperatures T below the mean field melting temper- 
ature T m by setting the parameter r ex T — T m . However, 
the free energy of Eq. |(3J) do not provide all the informa- 
tion on the system at finite temperatures and specially do 
not take into account strong thermal fluctuations which 
may lead to the melting of the hexagonal pattern. 

To go beyond mean field theory and take thermal fluc- 
tuations into account, we follow the usual procedure of in- 
troducing a coarse-grained effective Hamiltonian 11 where 
the energy for a particular configuration of the order pa- 
rameter is given by the corresponding mean field free- 
energy functional. Then instead of considering just the 
minimum energy this procedure allows for excitations ac- 
cording to a statistical weight provided by the free energy 
functional. In the present case we take as the effective 
Hamiltonian H[ip] = F[ip] with the corresponding parti- 
tion function 

Z= £ e- H ^ ksT , (3) 
MS)} 

where T is the temperature and the summation (func- 
tional integration) is taken over all configurations of tp(x). 



Physical quantities are defined in the usual way as ther- 
mal averages over iji configurations. The parameter r con- 
tained in H[ip] is taken to be a temperature-independent 
constant for the region below the mean field transition 
temperature T m in the absence of the pinning potential. 

We consider a pinning potential V(x, y) with square 
symmetry 

V(x,y) = Vo[cos(fc s a;) + cos(k s y)}, (4) 

where k s — 2Tr/a s is the wave vector of the pinning po- 
tential, which can represent, for example, and adsorbed 
layer on the (100) face of an fee crystal 1 ' 3 . We define 
the lattice misfit between the lattice system and pinning 
potential as S m = (fco — k s )/ko, where fco = 1 is the wave 
vector of the hexagonal periodic pattern^ of ijj(x,y) in 
the absence of the pinning potential, corresponding to a 
free lattice system. Thus, for S m 3> 0, the lattice sys- 
tem is under tensile strain while for <5 m <C it is under 
compression. 

III. MONTE CARLO SIMULATION 

The free energy associated with the partition function 
of Eq. ([3]), which now takes into account thermal fluc- 
tuations can not be obtained directly without approxi- 
mations. Therefore, to obtain the equilibration configu- 
rations and average quantities it is necessary to perform 
numerical simulations and for this purpose we used the 
MC method. Thermal averages over configurations of 
ij)(x,y) were performed using MC simulations on a dis- 
crete version of the effective Hamiltonian H [ip] . The field 
ij)(x,y) and pinning potential V(x,y) are defined on a 
square space grid, x = idx, y — jdy integers), with 
grid size L x L and grid spacing dx — dy and with pe- 
riodic boundary conditions. The Laplacians in Eq. @ 
were approximated by a simple discretization scheme 

d 2 t _ fi+hi + h-\,j - 2fi,j (fl 
dx* Ud ~ (Ax) 2 ■ lJ 

The standard Metropolis algorithm was used for simula- 
tions at a given temperature. At each grid site 
we attempt to change tpij by an small amount 
with probability min(l, e - AH / k sT^ usm g the Metropolis 
scheme, where AH is the resulting change in the configu- 
rational energy. Simulations were performed using, typ- 
ically, L ranging from 64 to 224, grid spacing dx = it/ 4 
and (4 — 2) x 10 6 MC passes for equilibration and equal 
number of passes for thermal averages. Because of the 
periodic boundary conditions, the length of the system 
L dx has to be a multiple of the lattice constant of the 
pinning potential, 2n/k s . This condition imposes a re- 
striction on the values of fc s when calculations are done 
for fixed L and dx. Thus, to be able to change k s con- 
tinuously near a commensurate phase while still keep- 
ing the same system size L, we choose a variable grid 
spacing dx — (fc°/fc s )(7r/4), where fcg is fixed such that 
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Ltt/A = n(27r/fc°) (n is an integer). The parameters r 
and /x were set to r = —1/4 and ju = —0.1875, corre- 
sponding to the hexagonal crystalline region in the orig- 
inal PFC model without pinning 8,9 . The temperature T 
is measured in units of Fqk'g 1 . 

Near the phase transitions where long equilibration 
times are required, we use the exchange MC method 
(parallel tempering) 13,14 . This method is known to re- 
duce significantly the critical slowing down near phase 
transitions. Recently, it has being used to study glassy 
incommensurate vortex lattices in 2D superconducting 
arrays 1 ^. In this method, many replicas of the system 
with different temperatures Tj in a range above and be- 
low the critical temperature are simulated simultane- 
ously and the corresponding configurations are allowed 
to be exchanged with a probability distribution satisfy- 
ing detailed balance. The exchange process allows the 
configurations of the system to explore the temperature 
space, being cooled down and warmed up, and the sys- 
tem can in principle escape more easily from metastablc 
minima at lower temperatures. Without the replica ex- 
change step, the method reduces to conventional MC 
simulations performed at different temperatures. The 
method was implemented by performing MC simulations 
as described above for each replica at different temper- 
atures, simultaneously and independently, for a few MC 
passes. Then exchange of pairs of replica configura- 
tions at temperatures Tj and Tj and energies Ei and 
Ej is attempted with probability min(l, exp(— A)), where 
A = (1/Tj — 1 /Tj )(Ej-Ei), using the Metropolis scheme. 
We typically used 10 6 MC passes for equilibration with 
up to 30 replicas and equal number of MC passes for 
calculations of average quantities. 

IV. RESULTS AND DISCUSSION 

A. Ground-state phase diagram 

The ordered structures in the ground state were ob- 
tained by simulations of the dissipative dynamical equa- 
tions for the phase field if}(x), as in the previous work 7 . 
The reason for using this method instead of the MC 
method described in Sec. Ill is that, at very low tem- 
peratures, MC simulations turned out to be computa- 
tionally less efficient due to the low acceptance rate for 
MC moves. The dissipative dynamics should evolve the 
system to the lower-energy state for arbitrary initial con- 
ditions. Therefore the determination of the final config- 
urations is equivalent to finding the ground state for the 
effective Hamiltonian in Eq. ([3]), since H\ip\ = F[ijj]. For 
the free-energy functional of Eq. ^ with non-conserved 
field the dynamical equations are 



The ground-state phase diagram was obtained by find- 
ing the final configurations where dtp /St = starting 




FIG. 1: Density plot of the phase field ip(x) showing commen- 
surate (C) and incommensurate (IC) structures in the ground 
state, depending on the amplitude of the pinning potential Vo 
and misfit parameter S m . (a) (1 x 1) C phase; (b) c(2 x 2) C 
phase; (c) (2 x 1) C phase; (d) hexagonal (full) IC phase and 
(e) IC phase with domain walls near the c(2 x 2) C phase. 

from different initial configurations, for different misfits 
5 m and amplitudes of the pinning potential Vq. Typi- 
cally, we start with a hexagonal structure which is the 
minimum energy state in the absence of the periodic po- 
tential (Vo = 0) and then increase slowly Vq to the its 
final value for fixed S m . The calculations are then re- 
peated for different initials conditions. The ground state 
is the final configuration corresponding to the lowest en- 
ergy. 

Three low-order commensurate structures were ob- 
tained as shown in Fig. [TJ a (1 x 1) commensurate (C) 
phase, where the peaks of ip(x) coincide with the pinning 
potential minima; a c(2 x 2) C phase, where they form a 
superstructure with periodicity twice that of the pinning 
potential along the principal directions and a c(2 x 1) 
phase where the superstructure has lattice periodicity 
twice the pinning potential along one of the directions. In 
addition, there are incommensurate (IC) phases, where 
ip(x) forms a hexagonal periodic structures incommen- 
surate with the pinning potential as in Fig. [ljd) or a 
structure of domain walls separating commensurate re- 
gions as in Fig. UJe). 

The results of extensive numerical calculations using 
the dynamical equation ([6]) are summarized in the phase 
diagram in Fig. [2] for —0.5 < S m < —0.225 and in Fig. 
[3] for —0.2 < 6 m < 0.40. The transitions between the 
different phases were determined from the change in the 
structure factor. 



B. Finite-temperature phase diagram 

We have studied the influence of thermal fluctuations 
and lattice mismatch near the simplest commensurate 
structures, the (1 x 1) and c(2 x 2) commensurate phases, 
using the MC method described in Sec. III. The phase 
diagrams were obtained by monitoring the behavior of 
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(8) 



FIG. 2: Ground state phase diagram in terms of the pinning 
strength Vo and mismatch 8 m for —0.5 < S m < —0.225 
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FIG. 3: Ground state phase diagram in terms of the pinning 
strength Vo and mismatch S m for —0.2 < 5 m < 0.4. 



the structure factor and specific heat as a function of 
temperature and lattice misfit. First, the structure factor 
S(k) was calculated from the positions Rj of the local 
peaks in the field ip(x,y) as 



S(k) 



N P 

<E 



N P 



-ig.(i? 3 — 



(7) 



where iV p is the number of peaks. To locate the peak po- 
sitions Rj for each configuration ofip(x, y) we have imple- 
mented a computer algorithm^ to find the local maxima 
in ip(x,y) based on a particle location algorithm used 
in digital image processing 17 . Alternatively, a structure 
factor can also be defined directly from the field ip(x,y) 



where ip(k) is the Fourier transform of ^{x). While the 
two expressions give similar results the former expression 
is better for characterizing a structural phase transition 
as it includes only the ordering of the lattice and does not 
simultaneously include fluctuations in the amplitude of 
the ip as the latter expression does. The results described 
here were obtained using the definition in Eq. (|7|). Sec- 
ond, the specific heat c was calculated directly from the 
average energy as c = (1/ L 2 )(d(H) /dT) and from the 
fluctuation relation 



1 



LP-T 2 



{(H 2 ) - (Hf ). 



(9) 



We checked that these two different ways of calculating 
c gave the same results, which indicates that the system 
is properly equilibrated at a given temperature. 

Figs. HJi and [4j) show the behavior of the scaled 
structure-factor peak S(k m )/N p and specific heat c as 
a function of temperature for the (1 x 1) and c(2 x 2) 
commensurate phases. Here k m is the magnitude of the 
wave- vector corresponding to the local maximum in S(k). 
For the (1 x 1) phase, S(k m )/N p in Fig. d£i decreases 
and broaden with temperature but does not disappear 
at the highest temperature where a liquid-like phase is 
expected. Also, there is no peak in the corresponding 
specific heat in Fig. [4)d. This indicates that the (lxl) 
phase does not melt into a liquid-like phase via a phase 
transition, instead there is a smooth crossover from a low 
temperature highly ordered phase to a high temperature 
phase where the pinning potential still induces some or- 
der in the peak pattern of ijj(x). This behavior is ex- 
pected on theoretical ground s 1 ' 18 since in this case the 
pinning potential has the same symmetry as the com- 
mensurate phase and acts as a constant external field on 
the displacement order parameter. On the other hand, 
for the c(2 x 2) commensurate phase, S(k m )/N p in Fig. 

drops sharply near the temperature where there is a 
peak in the corresponding specific heat in Fig. UJd and 
this behavior is associated with the melting of the ordered 
structure into a liquid-like phase. 

The phase diagrams obtained near the (lxl) and 
c(2 x 2) commensurate phases as a function of mis- 
fit, strength of pinning potential and temperature are 
shown in Figs. [5] and [6l respectively. In addition to 
the melting transitions, there are also commensurate- 
incommensurate transitions from the c(2 x 2) to IC phases 
and from c(2 x 2) to (1 x 1) phases which are identi- 
fied by the change in the peak patterns in the structure 
factor. A striking feature of the phase diagram is its 
topology, which is strong dependent on the type of C 
phase. This agrees qualitatively with theoretical predic- 
tions from simplified models^ 8 -. However, the present 
calculations were not sufficiently accurate to determine 
the joining of the transition lines near the incommen- 
surate phases in Figs. H>fa) and (b). Hence, we were 
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FIG. 4: Temperature dependence of the scaled structure- 
factor peak S(km)/N p (a) and specific heat c (b) for the (1 x 1) 
commnesurate phase ( 5m = 0, Vo = 0.10 ) and c(2 x 2) 
commensurate phase ( 8 m = —0.5, Vo = 0.275). Here k m is 
the wave vector of the corresponding ordered structure. 



FIG. 5: Phase diagrams near the (lxl) commensurate phase, 
(a) at fixed pinning strength Vo and (b) at fixed misfit <5 m . 



unable to investigate the theoretical prediction o 18 ' 19 for 
the existence of an intermediate liquid phase between the 
IC and C phases at finite temperature near the c(2 x 2) 
commensurate phase. 

In addition to the topology of the phase diagram, 
the nature of the phase transitions between the dif- 
ferent phases is of particular interest. We have in- 
vestigated in detail the critical behavior for the sim- 
plest case, corresponding to the melting transition of the 
c(2 x 2) commensurate phase as a function of the tem- 
perature. From symmetry considerations and Landau 
free-energy expansions 2 ^, one expects that the critical 
behavior should be in the Ising universality class. To 
verify if the phase field crystal model provides a correct 
description of this behavior, we have performed a finite- 
size scaling analysis of the melting of the c(2 x 2) phase 
as a function of temperature, for a fixed value of misfit 
<5 m = —0.5. Calculations were performed for increasing 
system sizes for the specific heat, structure factor and 
a suitable dimensionless Binder ratio^i. We define the 
Binder ratio Ul{T) as 

FIG. 6: Phase diagrams near the c(2 x 2) commensurate 
,, . . , 4 . phase, (a) at fixed pinning strength Vo and (b) at fixed misfit 

Ur=2- { \ P{km >\ I (10) Sm. 
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where the order parameter p{k m ) is the Fourier transform 
of the of density of local peaks in the field i/)(x, y) at Rj 



(a) 



j=i 



-ik m -Rj 



(ii) 



evaluated at the wave vector fc m corresponding to the lo- 
cal maximum of the structure factor S(k). The finite-size 
behavior of Ul provides an accurate determination of the 
critical temperature T c and an estimate of the thermal 
critical exponent v which characterizes the divergent cor- 
relation length, £ oc \T — T c \~ v , near the transitio n 21 ' 22 . 
In the high temperature disordered phase, the real and 
imaginary parts of p(k m ) fluctuate with a Gaussian dis- 
tribution leading to Ul — > while at low temperature 
there is long-range order with {p{k m )) ^ and therefore 
Ul — > 1 for L — > oo. At the critical temperature T c , the 
ratio U l becomes independent of the system size L and 
therefore plots of Ul (T) as a function of temperature for 
different systems sizes should cross at the same point, 
corresponding to the critical temperature T c of the sys- 
tem in the thermodynamic limit. In the scaling regime 
sufficiently close to T c , the dimensionless Ul(T) should 
satisfy the scaling form 



U L {T) = U{{T~T c )L 1 /») 



(12) 



where U (x) is a scaling function. Since the slope 
(dU L {T)/dT) evaluated at T c is proportional to L 1 ^, 
an estimate of v can be obtained 2 ^ from a log-log plot of 
this quantity against L. 

Fig. [7(a) shows the temperature behavior of Ul for dif- 
ferent system sizes. Although the curves do not cross pre- 
cisely at the same point, for the three largest system sizes 
the curves intersect approximately at T c w 0.03599(20). 
The lack of intersection at a common point should be 
due to statistical errors and corrections to scaling. Fig. 
[7b shows a log- log of SUl{T)/5T evaluated at the esti- 
mated T c against L from where we obtain \jv = 1.1(1). 
This estimate of v is in agreement with the exact value 
for the thermal exponent of the two-dimensional Ising 
model, v = 1. 

The finite-size behavior of the specific heat c is also 
consistent with the Ising universality class, where the 
specific heat exponent is a = 0. Fig. [5(a) shows the 
temperature behavior of the specific heat for the different 
systems sizes. From finite-size scaling, the maximum of c 
should scale as L 01 '", which corresponds to a logarithmic 
behavior C max oc logL, if a = 0. Fig. [5(b) shows that 
the linear-log plot of the specific heat maximum c rnax 
against L is indeed consistent with a logarithmic behav- 
ior for the three largest system sizes. 

The correlation function exponent 77 can also be esti- 
mated from the expected finite-size behavior of the struc- 
ture factor, which should scale as S(k m ) oc L 2 ~ v at T c . 
Fig. [9] shows a log-log plot of S(k m )/L 2 evaluated at 
the above estimated T c against L from where we find 
r) = 0.15. Unlike the above estimates of v and a, this es- 
timate of r\ is much smaller than the known exact value 
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FIG. 7: (a) Temperature dependence of the Binder ratio 
Ul (T) for different system sizes L, near the melting transition 
of the c(2 x 2) commensurate phase at 8 m = —0.5 and Vb = 
0.275; (b) Estimate of the thermal critical exponent v from 
the log- log plot of -J^Ul(T) at T c against L for the three 
largest system sizes. 



for the Ising model, i] = 0.25. This large discrepancy 
could be due to corrections to scaling. However, to take 
such corrections into account would require more accu- 
rate data and larger system sizes, which is beyond the 
scope of the present work. 



V. CONCLUSIONS 

In the present work we have studied the effects of ther- 
mal fluctuations in the phase field crystal model with an 
external pinning potential 7 using a non-conserved ver- 
sion of the model and Monte Carlo simulations. We have 
determined the phase diagram as a function of temper- 
ature and mismatch near low-order rational commensu- 
rate phases. The results show a rich phase diagram with 
commensurate, incommensurate and liquid-like phases 
with a topology dependent on the type of ordered struc- 
ture qualitatively consistent with predictions of simpli- 
fied models of pinned lattice systemsi^. In particu- 
lar, we find that the melting transition for the c(2 x 2) 
commensurate phase is consistent with the Ising uni- 
versality class, which is expected from analytical argu- 
ments based on symmetry considerations and Landau 
free-energy expansion 2 ^. Our results demonstrate that 
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the PFC model and the MC method employed here can 
be used to study specific lattice systems by adjusting 
the parameters of the model to match the experimen- 
tal structure factor— ^ and choosing a suitable pinning 
potential. 
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FIG. 8: (a) Temperature dependence of the specific heat c 
for different system sizes L, near the melting transition of the 
c(2 x 2) commensurate phase at <5 m = —0.5 and Vo = 0.275; 
(b) Specific heat maxima c max in a linear-log plot indicating 
a logarithmic behavior as a function of L for the three largest 
systems sizes. 
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